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We present a general analytic method to compute the number of metastable configurations as a 
function of the energy for a system of interacting Ising spins on the Bethe lattice. Our approach is 
based on the cavity method. We apply it to the case of ferromagnetic interactions, and also to the 
• binary and Gaussian spin glasses. Most of our results are obtained within the replica symmetric 

O ' ansatz, but we illustrate how replica symmetry breaking can be performed. 

PACS numbers; 



I. INTRODUCTION 



J3 ' Despite years of efforts, the nature of the glassy phase of finite dimensional spin glasses is still not clear. It is 
still debated whether the replica symmetry breaking (RSB) scheme proposed by Parisi to solve the mean field fully 
' , • connected Sherrington-Kirkpatrick model (SK), and implying the existence of an exponentially large number of pure 
^ [ states vifith an ultrametric structure, holds in some way for these systems. Part of the difficulty to settle the question 
Q i is the very poor analytical tractability of finite dimensional systems. 

^ ' Looking for more realistic ~ but still tractable - models than SK, much attention was recently paid to spin glasses 
on random graphs of finite connectivity. These models incorporate the short range nature of interactions, but without 
I ' the underlying geometry of the finite dimensional models. It has been shown that the replica symmetric (RS) solution 
^ of spin glasses on random graphs of finite connectivity is unstable at low temperature ||^, ^ . Unfortunately working 
out the RSB scheme is far more difficult than for the SK, because the glassy phase can no longer be characterized 
[ by the only two spin overlap {cTaUb) between two distinct replicas a and 6, but requires all of the multi spin overlaps. 
CO I Recently a new approach was suggested in based on the cavity method and population algorithms, that allows 
' for a numerical solution at the level of one step RSB (IRSB). The method can be virtually extended to any step of 
RSB, at the price of increasing computer resources. 

In this paper, we use similar ideas to shed a new light onto an old problem: the computation of the number of 
metastable configurations. A configuration is said to be metastable if its energy cannot be decreased by flipping a 

■ single spin. One expects, within the RSB scenario, that a consequence at zero temperature of the exponential number 
Ch \ of pure states, is an exponential number of metastable configurations. Beware however that these two concepts are not 

I ■ as obviously linked as it might seem: in systems with finite connectivity the zero temperature limit of the pure states 
■ O are not the metastable configurations, but are argued to be the configurations stable with respect to an arbitrary 
C , finite number of spin flips |§, 

Let us call configurational entropy Sc{E) the logarithm, divided by the number of spins, of the number of metastable 
configurations having an energy density equal to E. Lots of efforts have been devoted to the study of Sc for SK. 
^ More recently, attention turned to spin systems on random graphs. On random graphs with fixed finite connectivity, 
k>( ] annealed computations were carried out for the binary spin glass and the ferromagnet was addressed in |Q. On 

■ random graphs of fluctuating finite connectivity, the case of the ferromagnet was solved in |^ - also via a population 
d ' algorithm, but in a different context as ours ~, and the authors give hints how to address the quenched computation 

of spin glasses. 

Our method, which is quite general, enables us to recover all the above results. As new contributions, we carry out 
quenched computations in the case of the Gaussian or binary spin glass on random graphs of fixed finite connectivity, 
and we exemplify the implementation of IRSB. 

In the following we stick to the case of random graphs of fixed finite connectivity. For the sake of concision, we 
display no result about the fluctuating connectivity, but it is easy to generalize our method to this case. 

The layout of the paper is the following. In section || we give several definitions and we set the notations. In 



section III we derive our basic equations in the RS framework. In section [V we analytically solve these equations in 
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FIG. 1: Part of a Cayley tree 



the case of a ferromagnet. We compare our results with the microcanonical approach of 0|. In section |^ we turn to 
the binary spin glass, where the coupling constants are ±1. The population algorithm shows up there. In section VI 
we address the slightly more involved Gaussian spin glass. Eventually, in section VII we show how to perform IRSB, 
and illustrate the algorithm on the case of the spin glass with binary couplings. 



II. THE SYSTEM UNDER STUDY 

Following we call Bethe lattice a random graph with fixed connectivity equal to fc + 1, i.e. the number of edges 
incident to each vertex is exactly fc + 1. Locally such a graph has the structure of a Cayley tree, and can be obtained 
with the following procedure: starting from a root (or ancestor) '5, one builds a first generation of fc + 1 sons, and 
then successive generations of fc sons to each vertex, as displayed in Fig.(|l]). By contrast to the Cayley tree, the Bethe 
lattice has loops, but small ones are rare: the typical length of a loop is of order log A^. In practical case one uses 
the tree-like structure to write down recursive (or 'cavity' in the language of the physics of disordered systems [^) 
equations, and the existence of loops is enforced by saying that all the spins are equivalent, which is not true for the 
Cayley tree (in particular the Cayley tree has strong boundary effects). 

On each vertex A of the Bethe lattice stands a spin, and these spins interact through the Hamiltonian 

H — - ^ Ja,bo'ao'b, (1) 

{A,B} 

where the sum runs over all the edges. The coupling constants Ja,b may be fixed, like in the case of the ferromagnet, 
or they may be quenched random variables in the case of a spin glass. 

The local field acting on the spin A is Ha = J2b£V(A) Ja.bO'b, where the sum runs over the fc + 1 first neighbors of 
A. The metastable states are characterized by VA, Hacta > 0. Our aim is to compute the partition function restricted 
to the metastable configurations: 



Z = Tt 



\{Q{HA<yA) 



(2) 



where O is the Heaviside step- function such that 0(u) 
over all the possible values of all the spins. 



1 if u > 0, otherwise, and Tr stands for the summation 



III. THE CAVITY EQUATIONS 

Let us introduce, in close analogy with the construction of the infinite tree in pd| ], a simple labeling rule for the 
vertices of a tree: the sons of the root are labeled 1, 2, . . . fc. Then recursively the sons of a vertex labeled (ii . . . ip) 
are labeled (ii . . .ipl), . . . {ii . . Ap fc). If A = (zi . . . ip), we can define the norm of the symbol as |^| = p. 

The basic building block of the cavity solution is a branch of the Cayley tree, i.e. a subtree made of a given vertex 
different from the root and all its descendants. On Fig. (0) is drawn a branch rooted at Let us call it T^. 
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The key property of a branch is that it is structurally similar to each of its sub-branches made of a given vertex 
and all its descendants. This makes a recursive approach possible to study the thermodynamics of a branch. We will 
see afterwards how to use the results for the Bethe lattice. 



A. Computing the metastable states on a branch 

As it will be clear later, on a branch it is convenient to require that all the spins are stables, except the root. The 
Hamiltonian of a branch is 



7^$ 



JA.BO'ACrB- 



(3) 



The definition of the partition function is 



= Tr 



-fJHi. 



(4) 



By contrast to all the other spins in r$, which have fc + 1 neighbors, spin $ lacks one neighbor. In all the following 
the local field acting on such a spin with only k neighbors will be called a 'cavity' field and denoted by a small h (here 
We call p<t,{h, a) the joint density probability of and the value of spin <I>. One has 



p^{h,a) ^ ^Tr 



-0H. 



'S{a^ - cr) Jl Q{HAcrA)S{h^ - h) 



\A\>1 



(5) 



We call Ti, . . .Tfc the subtrees of T,^ engendered by points 1, . . . We can define their Hamiltonians similarly to 
Eq. (0), and write 



(6) 



Splitting the Tr, equation (^) can be restated as 

1 



4=1 



\A\>2 



p^{h,a) ^ ^el''^" Tr 5{h^ ~ h) Tr e-^^^ti TT e(iJ^a^) TT e(ff^a^). 

Z$ o-i,...crfc |-4|>2 

Now separate Hi into two contributions: 

k 



(7) 



(8) 



and enforce the fact that the first one is the cavity field acting on spin i in the absence of $ through the identities 

1= dht5(^^^Ji^tjaij - (9) 

Plugging ^ and into (|^) one gets 

k 

p^{h, a) = ^e'^'"^ Tr Sih^ ~ h)f] [ dh,e [{h, + J*,,a)a,] 

Z$ cri, ...CTfc ^J- J 



(10) 



\A\>2 



L4|>2 



i=l j 
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FIG. 2; A branch 




FIG. 3: Result of the merging (M2) of 2 branches 



Comparing to Eq. it turns out that the second hne in the above equation is nothing but the product of ZiPi{hi, ai) 
for i = 1, . . . k. So one eventuaUy obtains the following recursion relation: 

Piih,c7) = ^^-^-^^ Pha Sih^-h)Y\ [ dh,e[ih, + .U^,a)a,]p,{h,,a,). (11) 

Z<[. fl,...(Tfc J 



By averaging over the coupling constants and the random graphs (operation denoted by , one defines the proba- 
bility distributions V^'''\p) = S{pq,{h,a) —p). Equation (|ll| ) induces a functional relation between V^''^^ and the 

analogous probability distributions 'P^'^''^^ for z = 1, . . . /c. In the Bethe lattice all spins are required to be equivalent 
so we impose the condition 

Vi, V^^V^^V. (12) 

This yields a self-consistency equation for V. 



B. Performing measures on the Bethe lattice 

Having solved the thermodynamics of a branch, the idea of Q is to describe the Bethe lattice as the result of the 
merging of several branches. We define two merging procedures: 

(Ml) Take k + 1 branches, with roots labeled i = 1, . . .k and merge them onto a new spin The resulting tree is 
the one of Fig. |l|. 

(M2) Take 2 branches, with roots labeled 1 and 2, and merge them via a new bond of coupling constant J, without 
adding any site. See Fig. |[ 

Operation (Ml) is very similar to the iteration Eq. (|ll]): the differences are that k is to be substituted by -I- 1, 
and one must enforce the stability of the spin '5. After the merging, the joint density probability of the value of spin 
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^ and the local field acting on it is P^{H, a), given by: 



fe+i „ 

'./^"^eiHa) Tr 5{H^ - H)!] dhiO [{hi + J^,^a)a^] p^{hi, ai). (13) 

O-lv-O-fe + l ^-^ J 

1=1 



Summing over all the possible {H, a) and taking the logarithm one gets the variation of the free energy of the system 
during the merging: 



AFi = -iln 



V / dHef^"^e{Ha) Tr 5{H^ - H)f\ [ dhiO [{hi + J^,ia)ai] Pi{hi, at 

k+l ,. 

T> ef^"^''Q{H^a)W dhiQ[{hi + J<i.,ia)ai]pi{hi,ai) . 

<7,cri,...<7fe+i f J 



As far as operation (M2) is concerned, it is straightforward to derive the variation of the free energy 



AF. = --ln 



Tr e^-^"^"^ I dhxdh2@[{hx + J(j2)(Ti]Q[{h2 + Ja{)(j2]pi{hi,ar)p2{h2,(T2] 

Tl,0-2 J 



(14) 



(15) 



One can easily deduce the density of free energy F of the Bethe lattice from AFi and Ai^2- Assuming that each 
branch has a number of spins equal to N, the system after the merging (Ml) has N{k + 1) + 1 spins, and its free 
energy can by written in two ways: 



(A^(fc + 1) + 1)F = (fc + l)Fbranch + AFj , (16) 

where Fhranch is the free energy of one of the branches before the merging. For operation (M2), one has 

{2N)F = 2Fbranch+'^- (17) 

Elimination of Fbranch between the above equations yields 



F = AFi 



k + l 



AF2. 



(18) 



Note that the above derivation holds for any extensive and self-averaging observable. Thus the density of energy 
E of the Bethe lattice is 



E = AEi - 



k + l 



AEo 



(19) 



where AE'i, resp. AE2, is the variation of energy under the merging process (Ml), resp. (M2): 
AEi = -J2 [ dHP^{H,(j)H(j 

cr ■' 

Tr<,,ai......+i jy*<7e'5^*"e(^*'7)n^=iV dh, 6 [{h, + J^,,G)<j,]p,{h„ai) 



Tr„ 



and 



AE2 = -J 



a,a,,...a,+, e^"*''&{H^a)lliIl I dh^ S [{h^ + ,^a)a^] p^{h^ , a^) 

Tr^,,<,3 CTio-ae'^^"^'"^'"" / dhidh2 6 [{hi + Ja2)cri] 6 [(/i2 + Jcri)a2] Pi(/ii, cri)p2(/i2, ^2) 



(20) 



Tr,,,,^, e/3"''^i°'2 J dhidh2<d[{hi + Jg2)(Ti] Q[{h2 + J(Ji)<72]pi{hu(Ji)p2{h2.(y2) 
These results can be applied to a wide variety of models: up to now we have not specified the coupling constants. 



(21) 
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IV. THE FERROMAGNETIC CASE: AN EXACTLY SOLVABLE MODEL 



First, we concentrate on a simple case where all the above equations can be worked out analytically. The previous 
formalism allows for the computation of the metastable states of the k = 2 Ising ferromagnet. When the couplings all 
have the same absolute value \ J\ = 1, the cavity field can take only value in the set (—2, 0, 2). Thus integration with 
respect to the hi in Eq. ( pi] ) and other formulas reduces to a finite summation. Moreover, in the case of the Ising 
ferromagnet (i.e. J = 1) the system is homogeneous, therefore the cavity equations should not depend on the site index 
and we can define the joint probabilities po = p{h = —2, a = —1), pi = p{h = —2, a = 1), p2 = p{h = 0, cr = — 1), 
... , P5 = p{h = 2,cr = 1) of Eq. ( pi]) identically for all the lattice sites. A simple enumeration shows that Eq. (|Tl| ) 
reduces to the following system of equations: 



cpoe 

CP3 = 2piPi + 2p3P4 cpiC^"^^ 



p\ cpie '^f^ ^pI + 2pip3 + pI cp2 



2piP2 + 2piP4 



pl + 2pip3,+pl cp5e^'^ =pI, 



(22) 



together with the normalization condition X]i=oPi ~ ^- symmetries of the system fix some conditions on the 

values of pi. The relations po — Ps, Pi = Pi, P2 — Ps hold in the paramagnetic phase (PA). If the overall Z2 symmetry 
is spontaneously broken the system encounters a ferromagnetic (EM) phase characterized hy p2 = ps = P4 = P5 — 
(obviously there is also the solution with the opposite magnetization, for which po = pi = P2 = Ps = 0) ■ In the FM 
case the solution is given by 



Pa 



1 



1 



1 



p4/3 



Pi 



1 + e-4/3 



(23) 



The solution in the PA case is slightly more involved: 



Po 



c{x) 



c(x)D{xf 



D(x)4 D{xf D{xf^ 
64a;5 Sx^ Ax 



P2 



c{x)D{x)^ 
2^ 



(24) 



where x = (^) 3 and D{x) = \/x^^~~2x — x^. The problem of the iterative stabilit y o f both FM and PA solutions can 
be addressed studying the spectrum of the Jacobian of the system of equations ( p2| ) . It turns out that PA solution 

becomes unstable below a ^ = \og-^{2^) = 2.039091 and the FM one above Tf = = 2.885390. Once we 

have calculated the cavity fields we can measure the thermodynamic potential using the merging procedures explained 
in the former section. 

Let us consider first the FM solution. 



ASi = -3 



AS2 



AFi = -3 
2 



-1 Ai^2 



/3 



log(l 



-4/3 



Inserting Eq. ( ^ ) into Eq. (18) and Eq. (|l9| ) we readily obtain E = F = i. 
its ground state. 

The PA solution is worked analogously. 



(25) 



the system is completely frozen in 



AEi = - 



AEn = - 



Nijx) 
Di{x) 
N2{x) 

D2{x) 



AF,^--\og{D,{x)) 
AF2^~^log{D2{x)), 



(26) 



where 



Niix) 
D,{x) 
N2{x) 

D2{X) 



2p-i{?,eP + e'P) + (Spip2{ 



6K(e 
2pie 
2pie 



.3/3 + 2e'5) + 6^1^2(6^5 + 63^) 

6plp2{3e^^ + 2e^) + ^Pxpi{e^ + 36^'') + ^ple'^ 



2{pi+p2fe^ 



+ 2{pi+p2re' 



(27) 
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FIG. 4: Thermodynamic potentials E and F vs. T for the ferromagnet in the case k — 2 



and for the pi we use Eq. (p4|). In Fig. (Q) we display both energy and free energy as a function of the temperature. 
Starting in the low temperature phase, the system gets trapped in the ferromagnetic solution in which all spins are 
aligned. This solution is locally stable for the iteration Eq.(22|) up to temperature T/"*^. We define Tc = 2.228512 
as the temperature at which PA and FM solutions have the same energy. Exactly at this temperature we have the 
coexistence of the two phases of the system. Above Tc the PA solution acquires a lower free energy. In Fig. (^) we 
display the configurational entropy Sc{E) as a function of the energy E. The thin straight line is tangent to the 
curve exactly at Sc{Ec) where Ec = E{Tc) ~ —1.05613, while the dotted curve beneath the tangent is the result 
of the microcanonical computation of Dean et al. presented in . Note that the microcanonical branch can not be 
obtained in our canonical scheme. 



V. THE BINARY SPIN GLASS 



Let us now turn to a slightly more complicated case, the binary spin glass, where the coupling constants are 
quenched independent identically distributed (i.i.d.) random variables with the following law: 

\{J) = \{5{J-l) + 5{J+l)). (28) 

As in the above section, an important simplification occurs: the local field acting on a spin can have only a finite 
number of values. The cavity fields can have fc + 1 values: — /c, — fc + 2, . . . ,k, whereas the local field acting on a site with 
k + 1 first neighbors can have {k + 2) values: —k — 1, —k + 1, . . . fc + 1. But by contrast to the ferromagnet, the system 
is no longer homogeneous, and the self consistency equation (|l2|) is a tricky object. To tackle a similar equation, 
suggested a numerical solution via a population algorithm. The basic idea is that a probability distribution can be 
represented by a large collection of i.i.d. random variables distributed according to it. 

In our case, the algorithm works as follows: we use a large population of Af sites i = 1, . . . A/", each of which is 
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FIG. 5: Configurational entropy Sc vs. E for the ferromagnet in the case k = 2 (see text). 



characterized by a finite set of 2(fc + 1) numbers, the Pi{h, a). After random initiahzation, one iterates the following 
sequence: 

Algorithm I 

(i) Select fc + 2 sites ii, . . . ik+2 at random, and extract some couplings Ji, . . . Jk+i- 

(ii) Perform the merging (Ml) of the branches rooted at «i, . . .ifc+i onto the site 5* = that is compute AFi 
with Eq. ([^, AE'i with Eq. (|o|). 

(iii) Perform the merging (M2) of the branches rooted at ii and 12, that is compute IS.F2 with Eq. (|l5|), IS.E2 with 
Eq. (^l|), using for instance the couphng Ji. 

(iv) Update the population: perform the merging of the branches rooted at ii,...ik onto the site $ = ik+2, that is 
substitute pk+2{h, a) by the result of Eq. (pl|). 

The algorithm converges in a stochastic sense: after a sufficient number of iterations, the pi{h,a) are distributed 
according to V. So one can compute the average of the thermodynamic quantities with respect to the couplings by 
'time' averaging of the measures performed in (ii) and (iii). 

In the practical implementation of the algorithm we have scanned values of M ranging from 200 to 4000. A 
careful finite size scaling analysis shows that results are really mildly dependent on Af as soon as A/" > 1000, and 
the asymptotic extrapolation is always compatible, within statistical errors, with the biggest size we have simulated. 
The errors are calculated with standard binning procedure, discarding the first half of the simulation. In this way 
we have a complete control on the 'Hhermalization^^ of the algorithm. Most of the simulations have been performed 
also for the same system removing the stability condition on the local fields. We will always refer to these data as 
AC (all configurations) to distinguish them from the metastable (MS) set of data. For most of the simulations we 
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FIG. 6; Entropy vs. E for the binary spin glass in the case k = 2. Inset: zoom of the low-energy region. MC with, AC 
without one spin-flip stability condition. 



have performed from 10"^ x A/" up to 10* x J\f iterations of Algorithm I, which we have verified to be enough both for 
thermalization and numerical accuracy. 

The configurational entropy Sc is the MS entropy. 

One should also note that in the case of the Bethe lattice the choice of the distribution Eq. (|8|) is not crucial and the 
same results hold in the case of a distribution A( J) = p 5( J— 1) + (1— p)(5( J+1) regardless how small p is and even in the 
limiting case of a purely anti-ferromagnetic system. This observation allows for the analytic computation of the critical 
temperature and the critical energy of this model. Following the route specified in Sec. [V we calculated the equivalent 
of the system of equations ( p^ ) for purely anti-ferromagnetic couplings. The stability analysis of the paramagnetic 

solution gives the same temperature of the FM case, i.e. = log"^(2y^) and = E{Tc) = -15/14 « -1.071429. 

Below this temperature (energy) the replica symmetry is spontaneously broken. 

In figure (|^) we display the logarithm of the number of metastable states as a function of the energy. In the inset 
we zoom the results for the lowest energy. It is interesting to note that in the region E < —1.22 both AC and MS 
entropy seem to behave linearly on E, and a linear fit works really well. The two curves meet at E — —1.273(5) 
and S = 0.0172(5). This energy is clearly compatible with the replica symmetric value of the ground state Eq-^ = 
—23/18 = —1.277777 jl^. We have an excess of entropy at E^^ which is proportional to the density of spin fous 
in the ground state, i.e. the number of spins with zero local field - such a spin can be flipped without changing the 
energy. 

In figure (0) we compare our data with the annealed approximation presented in while in inset we zoom the 
difference Sann{E) — S{E) in the region E < E^. Above this value the annealed approximation is believed to be 
exact, and in fact our data fall on the analytic curve within the error bars, while below E^ the two curves split. It is 
interesting to note that the splitting, barely visible in the main panel for E > —1.2, seems to be exponential, at least 
near Ec (note the logarithmic scale on the y axis of inset). 

Let us stress here that it is straightforward to adapt the algorithm to the case of random graphs of fluctuating 
connectivity of finite mean c. In such a graph, given two points, there is an edge connecting them with probability 
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FIG. 7: Entropy vs. E for the binary spin glass in the case k — 2: comparison of the quenched and annealed results. Inset: 
zoom of the difference between the annealed and metastable (MS) entropy. 




c/N, and no edge with probability 1 — c/N . Thus the number of first neighbors of a given site is a random variable 
distributed according to a Poisson law of mean c. To implement this, one must extract the number fc + 1 according 
to this law in step (i). The rest is unchanged, except the important fact that step (iv) is no longer to be performed 
with k branches, but with k + 1 branches, like (ii) and that the fc + 1 term in Eqs. (Es) and (19) should be replaced 
byc§. 



VI. THE GAUSSIAN SPIN GLASS 



Our method can tackle even more complicated cases, such as the Gaussian spin glass, in which the coupling constants 
are quenched i.i.d random variables whose law is a Gaussian with unit variance and zero mean. 

The local fields are now continuous variables. We choose to represent the probability distribution Pi{h, a) on a site 
by a large population of C couples (/i, a) distributed according to Pi. 

So on each site i — l,...Af, we have a population (/ij', erf), — !,...£. The algorithm's layout is similar to 
algorithm I: 

Algorithm II 

(i) Select fc + 2 sites ii, . . . ik+2 a-t random, and extract some couplings Ji, . . . Jk+i- 

. . ik+i onto the site 5* = ik+2- In the context of this algorithm. 



(ii) Perform the merging of the branches rooted at ii , 
the Eq.(|lJ) and (|2^) become 



■In 



fe+i 



e^^>e(i?^a)n0[(/^r, + vK 



U=l CT=±1 



(29) 
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AEi 



(30) 



where Jj^ij 



(iii) Perform the merging of the branches rooted at «i and 12, using for instance the coupUng Ji. The formulas ( p^ ) 
and (I21I) read: 



E 



(31) 



AE-2 



EL 



--e [{hi + Jiai)ai] e [{hi + JiO<^: 



2. 



.=1 eP'J^<^n-'^Q[{hl+J,al)al] 6 [(/j- 



(32) 



(iv) Update the population, performing the merging of the branches rooted at ii, . . .ik onto the site $ = ik+2' one 
wants to substitute the (ft,^, trj), j/ = 1, . . . £ with a new population distributed according to pi{h, a) of Eq. ([ll]). 
This is performed in two steps 

— First one builds a list {h'^ ,a'^),iy — 1,...£ as follows. Start with = 1. Enter a loop in which you 
extract randomly an ly' G {!,...£} and a spin value a — ±1, until you have the stability condition 
Vj S {1, ■ • ■ k}, {hi + Jja)<7l > 0. Then set h" = J2j=i k ^j'^ij — ^- Increment v, and enter the 
loop again. 

— The list {\f , a'^) is not the one to overwrite the (/i$, cr^), because one must enforce the presence of the factor 
g/3/i(T jj-^ (pT|). Thus the elements of the list are to be reweighted: a suitable new population {h'^,a^),iy = 
1, ... £ is obtained by repeating C times the process of picking an element in the list {h'^, a'^),v = 1, . . .£ 
with a probability proportional to e^'^' . 

We implemented the above algorithm in the case k — 2, using the values Af ~ 1000 and C = 1000. The results are 
presented in Fig. (||). By contrast with the binary spin glass, the ground state is not degenerated. So both entropies 
AC and MS go to at the energy density of the ground state Egs- Using the algorithm of |^ at low temperature 
and extrapolating the results to T = 0, one finds Egs = —1-12(1), which is compatible with our data. 

VII. BREAKING THE REPLICA SYMMETRY 

Let us now see how our method can be generalized to perform RSB. 



A. General considerations 



So far we have assumed there is a single pure state, stable under the iteration process. It is well known that this is 
actually not true for the spin glass on the Bethe lattice |^ . One should allow for an infinity of pure states as described 
by the hierarchical continuous RSB scheme proposed by Parisi [lO|] . As a first approximation, one can implement the 
IRSB scheme. Briefly speaking, one assumes that there exists an infinity of pure states labeled by a = 1, • • • + 00. The 
free energies of the states on one branch are independent identically distributed random variables, with an exponential 
density 



p{F) = exp{f3x{F - F^)) (33) 

where F^ is a reference free energy, and x is Parisi's parameter [|l^. In this approach one computes the average 
values of the observables as a function of x. The physical value of x is the one which maximizes the free energy. 
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FIG. 8: Entropy vs. E for the Gaussian spin glass. Inset: Energy vs. T = 1//3. Lines are power-law fits of the low-T region. 
Note that both curves extrapolate at T = to the same value (see text). 



If one considers a branch rooted at site in each pure state a the probability distribution p%{h, a) is different, and 
the iteration Eq. ( pi] ) is valid only inside a given pure state. The self-consistency equation to be solved is consequently 
more complicated than in the replica symmetric case. Let us write Pi(/i, ct) = (pf (^, tT))^^!..... One is interested in 

the set of probability distributions Q^''''^^(p) ~ (5(p$(/i, ti 



p), which are functionals of the analogous quantities 

Q^^''^\ and one asks for all of them to be equal. A comprehensive theoretical description of the extension of the cavity 
method to IRSB can be found in Q. 

Let us turn to the description of the algorithm. We use a population of M sites i — 1, . . .A/", M states a — 1, ... M. 
On each site, each state is characterized by its own distribution pf{h, a). Depending on the distribution of the coupling 
constants, pf{h, a) is represented as a set of 2(fc + 1) numbers (case of binary couplings), or again as a population of 
TZ couples (ft, cr) (case of Gaussian couphngs). 

The RS algorithm I is nested into the present algorithm: it is used to implement the iteration ( pi]) and compute 
the expectation values of the observables inside each state. Then a meta-algorithm takes into account all the states, 
with appropriate weights, to update the population and compute the global expectation values. Two new observables 
are required with respect to the RS case. First the variation of the free energy during the iteration process described 
by Eq. (O): 



AF.te 
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In 



Tr ^Phita 
(T,ai,...ah 



J dhiQ [{h, + J^sa)ai] Pi{hi, a^) 

i—l 

Second the derivative of the free energy F with respect to x. By derivation of Eq. ( p^ and (|l5|), one gets 

dx X 2 ' 



(34) 



(35) 
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where 



(36) 



and similarly for (P. 

An iteration of the algorithm goes as follows: 

Algorithm III 

(i) Perform step (i) of algorithm I. 

(ii) For each state a, perform steps (ii), (iii) and (iv) (computing en passant AF^^^.), of algorithm I. One gets 
quantities bearing a as a superscript: AF" , .... 



The states with low AF^^^ have to be favored. Thus 



(iii) Reweight the states: this is a crucial step, motivated in 
one picks up M times an element in the list of distributions p%,a = 1,...M with a probability proportional 
to exp{—f3xAF°^^^). The resulting list overwrites the one at site $. (Note that the elements one picks up are 
composite objects, i.e. either a set of 2(fc + 1) numbers or a population of TZ couples). 

(iv) Compute the global average values of the observables AFi .... The formulas can be found in 0] : 



(3x 



In 



M ^ 



(37) 



Ai^i = . (38) 

and their obvious analogues for AF2 and AE2. 

Note that it is essential that the sites and couplings be the same for all the a in step (ii) above. 

In practice, we found that the implementation of this algorithm in the case of the Gaussian couplings requires too 
many computer resources to reach a satisfying accuracy. So we limited ourselves to the case of binary couplings. 

B. Application to the binary spin glass 

We considered the binary spin glass of section at temperature T — 0.5. Our implementation used the values 
M — 1000 and A4 — 1000. We computed dF/dx for several values of x. The result is presented is Fig. (^. To 
determine the value x* of x where dF/dx is zero, we fitted the curve by a polynomial of degree 4, and found its 
roots. We got x* = 0.19(1). Then we measured the thermodynamic quantities for x = x* ^ taking into account the 
uncertainty on x* : 

F = -1.280 ±0.001, £;= -1.265 ±0.001, S* = 0.029 ± 0.001. (39) 
This is to be compared to the output of the replica symmetric algorithm of section ^: 

F= -1.2816 ±0.0008, S = -1.2685 ± 0.0006, 5* = 0.0274 ± 0.0008. (40) 

The overall improvement is small, particularly on F. It is more obvious on E and S but still only of the order of 
10-3. 
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FIG. 9: dF/dx vs. x for the binary spin glass in the case k = 2 aX T=0.5; the continuous line is the polynomial fit of degree 4 



VIII. CONCLUSION 

In this paper we have described a general procedure to compute the number of 1-spin flip stable configurations 
on the Bethe lattice. Given some integer fc > 2, the method - at least conceptually - can be easily generalized 
to compute the number of fc-spin fiips stable configurations, i.e. whose energies can not be decreased by flipping a 
number of spins ranging from 1 to k. The practical difficulty is that the recursion relations can no longer involve 
only quantities related to the root of a branch, but must also take into account the k — 1 successive generations of 
spins. In the case fc = 2 it remains quite straightforward to write down the equations, but we have not worked out 
their solution. It might be interesting to do so in order to clarify the nature of the zero temperature limit of the pure 
states: recently Mezard and Parisi Q presented a computation at IRSB level of the number of locally ground states 
(LGS), i.e. configurations stable with respect to fc-spin flips with the number fc going to infinity with the size of the 
system in some unprecised way, and they found some surprising features. There is a need for more precise definitions, 
and one of the points at issue is how these LGS are related to fc-spin flips stable configurations when fc -\-oo. 
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